#### installing packages if not installed ####

list.of.packages <- c('dplyr',
                      'estimatr',
                      'tidyverse',
                      'haven',
                      'readxl',
                      'ggthemes',
                      'rdrobust',
                      'gridExtra',
                      'lubridate')
new.packages <- 
  list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

#### purpose: analyze terry stop data #### 

#### attach packages ####

suppressPackageStartupMessages(
  
  {
    
    library(dplyr)
    library(estimatr)
    library(tidyverse)
    library(haven)
    library(readxl)
    library(ggthemes)
    library(rdrobust)
    library(gridExtra)
    library(lubridate)
    
  }
  
)

#### load data #### 

# loading stop data 

load("data/terry.RData")

terry = terry %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                    call_type == "SCHEDULED EVENT (RECURRING)" | 
                      call_type == "ALARM CALL (NOT POLICE ALARM)", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE", 1, 0))

# loading 911 call data
 
load(file = "data/calldata.RData")

call = 
  call %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)" | 
                                 call_type == "PROACTIVE (OFFICER INITIATED)" | 
                                 call_type == "POLICE (VARDA) ALARM", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE" | 
                    call_type == "HISTORY CALL (RETRO)" | 
                    call_type == "IN PERSON COMPLAINT", 1, 0))

# gen time data

call$time_arrived2 = parse_date_time(call$time_arrived, "%m/%d/%Y %H:%M:%S %p")
call$time_queued2 = parse_date_time(call$time_queued, "%m/%d/%Y %H:%M:%S %p")

call$time_arrived2 = 
  as.POSIXct(call$time_arrived, format="%m/%d/%Y %H:%M:%S %p", tz="UTC")
call$time_queued2 = 
  as.POSIXct(call$time_queued, format="%m/%d/%Y %H:%M:%S %p", tz="UTC")
call$response_time = (call$time_arrived2 - call$time_queued2) / 60
the_index = call$response_time < 0
call = call %>% filter(!the_index)

# loading in crime data 

load(file = "data/crimedata.RData")

crime = 
  crime %>% 
  mutate(date = as.Date(paste(substring(report_date, 7, 10),
                              substring(report_date, 1, 2),
                              substring(report_date, 4, 5), sep = "-"))) %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(cat_society = ifelse(crime_against_category == "SOCIETY", 1, 0),
         cat_property = ifelse(crime_against_category == "PROPERTY", 1, 0),
         cat_person = ifelse(crime_against_category == "PERSON", 1, 0),
         cat_notcrime = ifelse(crime_against_category == "NOT_A_CRIME", 1, 0)) %>% 
  filter(offense != "Identity Theft")

#### descriptive plot #### 

# outcomes of interest are: 
# 1) stops
# 2) officer initiated 911 calls 
# 3) hit rates
# 4) arrest rates
# 5) response times
# 6) against person
# 7) against property
# 8) against society 
# 4x3

# a) stops

terry_desc = terry %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm)) %>% 
  filter(date >= as.Date("2020-03-23"))
terry_desc$blmrv = 
  terry_desc$date - min(terry_desc$date[terry_desc$blm == 1])
terry_desc = terry_desc %>% filter(blmrv <= 67)

descp1 = terry_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Terry Stops",
       title = "A. Terry Stops") + 
  theme_tufte(base_size = 9)

# b) officer initiated calls 

call_desc = call %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count),
            blm = mean(blm)) %>% 
  filter(date >= as.Date("2020-03-23"))
call_desc$blmrv = 
  call_desc$date - min(call_desc$date[call_desc$blm == 1])
call_desc = call_desc %>% filter(blmrv <= 67)

descp2 = call_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .3,
             alpha= .3) + 
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Officer-Intiated Calls",
       title = "B. Officer-Initiated Calls") + 
  theme_tufte(base_size = 9)

# frisk rate

terry_hit_desc = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(hit = ifelse(stop_resolution == "Arrest" | 
                        stop_resolution == "Citation / Infraction" | 
                        stop_resolution == "Offense Report" | 
                        stop_resolution == "Referred for Prosecution", 1, 0)) %>% 
  summarize(hit = mean(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            frisk = mean(frisk, na.rm = TRUE)) %>%
  filter(date >= as.Date("2020-03-23"))

terry_hit_desc$blmrv = 
  terry_hit_desc$date - min(terry_hit_desc$date[terry_hit_desc$blm == 1])
terry_hit_desc = terry_hit_desc %>% filter(blmrv <= 67)

descp3 = terry_hit_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = frisk),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = frisk, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Frisk Rates",
       title = "Frisk Rates") + 
  theme_tufte(base_size = 9)

# c) hit rates

terry_hit_desc = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(hit = ifelse(stop_resolution == "Arrest" | 
                        stop_resolution == "Citation / Infraction" | 
                        stop_resolution == "Offense Report" | 
                        stop_resolution == "Referred for Prosecution", 1, 0)) %>% 
  summarize(hit = mean(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>%
  filter(date >= as.Date("2020-03-23"))

terry_hit_desc$blmrv = 
  terry_hit_desc$date - min(terry_hit_desc$date[terry_hit_desc$blm == 1])
terry_hit_desc = terry_hit_desc %>% filter(blmrv <= 67)

descp3 = terry_hit_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = hit),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = hit, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Hit Rates",
       title = "C. Hit Rates") + 
  theme_tufte(base_size = 9)

# d) arrest rates

terry_arr_desc = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(arrest = ifelse(stop_resolution == "Arrest", 1, 0)) %>% 
  summarize(arrest = mean(arrest, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>%
  filter(date >= as.Date("2020-03-23"))

terry_arr_desc$blmrv = 
  terry_arr_desc$date - min(terry_arr_desc$date[terry_arr_desc$blm == 1])
terry_arr_desc = terry_arr_desc %>% filter(blmrv <= 67)

descp4 = terry_arr_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = arrest),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = arrest, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Arrest Rates",
       title = "D. Arrest Rates") + 
  theme_tufte(base_size = 9)

# e) response times

call_rt_desc = call %>% 
  filter(initiate_civ == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(response_time = mean(response_time),
            blm = mean(blm)) %>% 
  filter(date >= as.Date("2020-03-23"))

call_rt_desc$blmrv = 
  call_rt_desc$date - min(call_rt_desc$date[call_rt_desc$blm == 1])
call_rt_desc = call_rt_desc %>% filter(blmrv <= 67)

descp5 = call_rt_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = response_time),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = response_time, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Response Time (in Minutes)",
       title = "E. Response Times") + 
  theme_tufte(base_size = 9)

# f) black/white disparity

terry_bw_desc = 
  terry %>% 
  mutate(count = 1,
         r_unk = ifelse(subject_race == "Unknown" | subject_race == "-", 1, 0),
         r_wht = ifelse(subject_race == "White", 1, 0),
         r_asn = ifelse(subject_race == "Asian", 1, 0),
         r_blk = ifelse(subject_race == "Black or African American", 1, 0)) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE),
            r_wht = sum(r_wht, na.rm = TRUE),
            r_blk = sum(r_blk, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  mutate(wht_rate = (r_wht / (.673 * 737015)) * 10000,
         blk_rate = (r_blk / (.073 * 737015)) * 10000) %>% 
  mutate(rr_bw = blk_rate / wht_rate) %>% 
  mutate(rr_bw = ifelse(rr_bw == Inf | is.na(rr_bw), NA, rr_bw)) %>% 
  filter(date >= as.Date("2020-03-23"))

terry_bw_desc$blmrv = 
  terry_bw_desc$date - min(terry_bw_desc$date[terry_bw_desc$blm == 1])
terry_bw_desc = terry_bw_desc %>% filter(blmrv <= 67)

descp6 = terry_bw_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = rr_bw),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = rr_bw, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Black/White Stop Rate Ratio",
       title = "F. Black/White Stop Rate Ratio") + 
  theme_tufte(base_size = 9)

# f) against person, property, society

crime_desc = crime %>% 
  group_by(date) %>% 
  summarize(cat_person = sum(cat_person),
            cat_property = sum(cat_property),
            cat_society = sum(cat_society), 
            blm = mean(blm)) %>% 
  filter(date >= as.Date("2020-03-23"))

crime_desc$blmrv = 
  crime_desc$date - min(crime_desc$date[crime_desc$blm == 1])
crime_desc = crime_desc %>% filter(blmrv <= 67)

descp7 = crime_desc %>% 
  filter(date >= as.Date("2014-01-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = cat_person),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = cat_person, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Crimes (Against Person)",
       title = "G. Crime (Against Person)") + 
  theme_tufte(base_size = 9)

descp8 = crime_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = cat_property),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = cat_property, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Crimes (Against Property)",
       title = "H. Crime (Against Property)") + 
  theme_tufte(base_size = 9)

descp9 = crime_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = cat_society),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = cat_society, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Crimes (Against Society)",
       title = "I. Crime (Against Society)") + 
  theme_tufte(base_size = 9)

descp_grob = arrangeGrob(descp1, descp2, descp3,
            descp4, descp5, descp6,
            descp7, descp8, descp9, ncol = 3)

ggsave(plot = descp_grob, filename = "pics/descfin.png", width = 8, height = 6.25)

#### descriptive plot (removing CHAZ) #### 


terry_desc = terry %>% 
  filter(initiate_off == 1) %>% 
  filter(beat != "E1") %>% 
  filter(beat != "E2") %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm)) %>% 
  filter(date >= as.Date("2020-03-23"))
terry_desc$blmrv = 
  terry_desc$date - min(terry_desc$date[terry_desc$blm == 1])
terry_desc = terry_desc %>% filter(blmrv <= 67)

descp1 = terry_desc %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Terry Stops",
       title = "A. Terry Stops") + 
  theme_tufte(base_size = 9)

# hit rates


terry_hit_desc = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(hit = ifelse(stop_resolution == "Arrest" | 
                        stop_resolution == "Citation / Infraction" | 
                        stop_resolution == "Offense Report" | 
                        stop_resolution == "Referred for Prosecution", 1, 0)) %>% 
  summarize(hit = mean(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

terry_hit_desc$blmrv = 
  terry_hit_desc$date - min(terry_hit_desc$date[terry_hit_desc$blm == 1])

descp3 = terry_hit_desc %>%
  ggplot() + 
  geom_point(aes(x = date, y = hit),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = hit, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Hit Rates",
       title = "C. Hit Rates") + 
  theme_tufte(base_size = 9)


# over time plot, effectiveness, assessing the raw data; makes Appendix figure D6
terry_arr_ts = 
  terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = sum(arrest, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01"))


terry_arr_ts$blmrv = 
  terry_arr_ts$date - min(terry_arr_ts$date[terry_arr_ts$blm == 1])

rdout = 
  with(terry_arr_ts, rdrobust(x = blmrv, y = arrest, p = 1, kernel = "tri", all = TRUE))

# hits

terry_hit_ts = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(hit = ifelse(stop_resolution == "Arrest" | 
                        stop_resolution == "Citation / Infraction" | 
                        stop_resolution == "Offense Report" | 
                        stop_resolution == "Referred for Prosecution", 1, 0)) %>% 
  summarize(hit = sum(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>%
  filter(date >= as.Date("2019-06-01"))

terry_hit_ts$blmrv = 
  terry_hit_ts$date - min(terry_hit_ts$date[terry_hit_ts$blm == 1])

rdout2 = 
  with(terry_hit_ts, rdrobust(x = blmrv, y = hit, p = 1, kernel = "tri", all = TRUE))

arrest_num = 
  terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = sum(arrest, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = arrest),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = arrest, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "A. Arrests") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout$coef[1], 2),
                                  "\nSE: ",
                                  round(rdout$se[3], 2),
                                  "\np < 0.05",
                                  "\npre-BLM SD: ",
                                  round(sd(terry_arr_ts$arrest[terry_arr_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = 6,
           family = "serif",
           size = 2.25) + 
  theme_tufte()

hit_num = terry_hit_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = hit),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = hit, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "B. Hits") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout2$coef[1], 2),
                                  "\nSE: ",
                                  round(rdout2$se[3], 2),
                                  "\np < 0.001",
                                  "\npre-BLM SD: ",
                                  round(sd(terry_hit_ts$hit[terry_hit_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = 6,
           family = "serif",
           size = 2.25) + 
  theme_tufte()

grobEff = arrangeGrob(arrest_num, hit_num, ncol = 2)
ggsave(plot = grobEff, filename = "pics/eff2.png", width = 8, height = 2.5)

#### rdit: cct ####

# polynomial: 0, 1, 2, 3
# kernel: 0, 1, 2, 3

# terry stop dataset

terry_cct = terry %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm)) %>% 
  mutate(count_stops = count)

terry_cct$blmrv = 
  terry_cct$date - min(terry_cct$date[terry_cct$blm == 1])

# call dataset

call_cct = call %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count),
            blm = mean(blm)) %>% 
  mutate(count_calls = count)

call_cct$blmrv = 
  call_cct$date - min(call_cct$date[call_cct$blm == 1])

# hit rates

terry_hit_cct = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(hit = ifelse(stop_resolution == "Arrest" | 
                        stop_resolution == "Citation / Infraction" | 
                        stop_resolution == "Offense Report" | 
                        stop_resolution == "Referred for Prosecution", 1, 0)) %>% 
  summarize(hit = mean(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

terry_hit_cct$blmrv = 
  terry_hit_cct$date - min(terry_hit_cct$date[terry_hit_cct$blm == 1])

# arrest rates

terry_arr_cct = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(arrest = ifelse(stop_resolution == "Arrest", 1, 0)) %>% 
  summarize(arrest = mean(arrest, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

terry_arr_cct$blmrv = 
  terry_arr_cct$date - min(terry_arr_cct$date[terry_arr_cct$blm == 1])

# response times

call_rt_cct = call %>% 
  filter(initiate_civ == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(response_time = mean(response_time),
            blm = mean(blm)) %>% 
  filter(date >= as.Date("2020-03-23"))

call_rt_cct$blmrv = 
  call_rt_cct$date - min(call_rt_cct$date[call_rt_cct$blm == 1])

# black/white disparity

terry_bw_cct = 
  terry %>% 
  mutate(count = 1,
         r_unk = ifelse(subject_race == "Unknown" | subject_race == "-", 1, 0),
         r_wht = ifelse(subject_race == "White", 1, 0),
         r_asn = ifelse(subject_race == "Asian", 1, 0),
         r_blk = ifelse(subject_race == "Black or African American", 1, 0)) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE),
            r_wht = sum(r_wht, na.rm = TRUE),
            r_blk = sum(r_blk, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  mutate(wht_rate = (r_wht / (.673 * 737015)) * 10000,
         blk_rate = (r_blk / (.073 * 737015)) * 10000) %>% 
  mutate(rr_bw = blk_rate / wht_rate) %>% 
  mutate(rr_bw = ifelse(rr_bw == Inf | is.na(rr_bw), NA, rr_bw)) 

terry_bw_cct$blmrv = 
  terry_bw_cct$date - min(terry_bw_cct$date[terry_bw_cct$blm == 1])

# against person, property, society

crime_cct = crime %>% 
  group_by(date) %>% 
  summarize(cat_person = sum(cat_person),
            cat_property = sum(cat_property),
            cat_society = sum(cat_society), 
            blm = mean(blm))

crime_cct$blmrv = 
  crime_cct$date - min(crime_cct$date[crime_cct$blm == 1])

# loop 

polys = c(0, 1, 2, 3)
kernz = c("tri", "uni", "epa")
datasets = list(terry_cct %>% as.data.frame,
                call_cct %>% as.data.frame,
                terry_hit_cct %>% as.data.frame,
                terry_arr_cct %>% as.data.frame, 
                call_rt_cct %>% as.data.frame,
                terry_bw_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame)
outcomes = 
  c("count_stops", "count_calls", "hit", 'arrest', 'response_time',
    "rr_bw", 'cat_person', "cat_property", "cat_society")

dataset_list = as.list(rep(NA, length(datasets)))

for (i in 1:length(datasets)) {
  
  print(paste0("I iteration ", i))
  polys_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J iteration ", j))
    
    kerns_matrix = matrix(NA, ncol = 6, nrow = length(kernz))
    
    for (k in 1:length(kernz)) {
      
      print(paste0("K iteration ", k))
      out = rdrobust(y = as.numeric(datasets[[i]][, outcomes[i]]),
               x = datasets[[i]]$blmrv, 
               p = polys[j],
               kernel = kernz[k])
      
      kerns_matrix[k, 1] = out$coef[3]
      kerns_matrix[k, 2] = out$se[3]
      kerns_matrix[k, 3] = out$pv[3]
      kerns_matrix[k, 4] = outcomes[i]
      kerns_matrix[k, 5] = polys[j]
      kerns_matrix[k, 6] = kernz[k]
      
    }
    
    polys_list[[j]] = 
      kerns_matrix %>% 
      as.data.frame %>% 
      `colnames<-` (c("est", "se", "pv", "outcomes", "poly", "kern")) %>% 
      mutate(est = as.numeric(as.character(est)),
             se = as.numeric(as.character(se)),
             pv = as.numeric(as.character(pv)),
             outcomes = as.character(outcomes),
             poly = as.numeric(as.character(poly)),
             kern = as.character(kern))
    
  }
  
  dataset_list[[i]] = polys_list
  
}

cct_estimates = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(kernpoly = paste0(kern, ", ", poly)) %>% 
  mutate(kernpoly = factor(kernpoly, 
                           levels = c("tri, 0", "tri, 1", "tri, 2", "tri, 3",
                                      "uni, 0", "uni, 1", "uni, 2", "uni, 3",
                                      "epa, 0", "epa, 1", "epa, 2", "epa, 3"),
                           labels = c("Tri.,\n0", "Tri.,\n1", "Tri.,\n2", "Tri.,\n3",
                                      "Uni.,\n0", "Uni.,\n1", "Uni.,\n2", "Uni.,\n3",
                                      "Epa.,\n0", "Epa.,\n1", "Epa.,\n2", "Epa.,\n3"))) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_stops", "count_calls", "hit", "arrest",
                                      "response_time", "rr_bw", "cat_person",
                                      "cat_property", "cat_society"),
                           labels = c("A. Terry Stops", "B. Officer-Initiated Calls", "C. Hit Rates",
                                      "D. Arrest Rates", "E. Response Time", "E. B/W Rate Ratio",
                                      "F. Crime (Against Person)", "G. Crime (Property)", "H. Crime (Against Society)"))) 

cct_estimates = 
  cct_estimates %>% filter(outcomes != "E. Response Time")

cct_estimates %>% 
  ggplot() + 
  geom_point(aes(x = kernpoly, y = est)) +
  geom_errorbar(aes(x = kernpoly,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                size = .2,
                width = 0) + 
  geom_errorbar(aes(x = kernpoly,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                size = .4,
                width = 0) + 
  geom_hline(yintercept = 0) + 
  facet_wrap(~outcomes, scales = "free", ncol = 4) + 
  labs(x = "Kernel, Polynomial",
       y = "RDiT Coefficient (BLM)") + 
  theme_tufte(base_size = 9) + 
  theme(axis.text.x = element_text(size = 3.5))

ggsave(plot = last_plot(), filename = "pics/cctres.png", width = 8, height = 3.5)

#### rdit: 10-65 day bandwidth ####

# loop 

polys = c(0, 1, 2, 3)
kernz = c("tri", "uni", "epa")
datasets = list(terry_cct %>% as.data.frame,
                call_cct %>% as.data.frame,
                terry_hit_cct %>% as.data.frame,
                terry_arr_cct %>% as.data.frame, 
                call_rt_cct %>% as.data.frame,
                terry_bw_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame)
outcomes = 
  c("count_stops", "count_calls", "hit", 'arrest', 'response_time',
    "rr_bw", 'cat_person', "cat_property", "cat_society")

dataset_list = as.list(rep(NA, length(datasets)))

for (i in 1:length(datasets)) {
  
  print(paste0("I iteration ", i))
  polys_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J iteration ", j))
    
    # kerns_matrix = matrix(NA, ncol = 6, nrow = length(kernz))
    kerns_list = as.list(rep(NA, length(kernz)))
    
    for (k in 1:length(kernz)) {
      
      print(paste0("K iteration ", k))
      
      bw_matrix = matrix(NA, ncol = 7, nrow = length(10:100))
      
      for (l in 1:length(10:100)) {
        
        print(paste0("L iteration ", l))
        out = rdrobust(y = as.numeric(datasets[[i]][, outcomes[i]]),
                       x = datasets[[i]]$blmrv, 
                       p = polys[j],
                       kernel = kernz[k],
                       h = (10:100)[l])
        
        bw_matrix[l, 1] = out$coef[3]
        bw_matrix[l, 2] = out$se[3]
        bw_matrix[l, 3] = out$pv[3]
        bw_matrix[l, 4] = outcomes[i]
        bw_matrix[l, 5] = polys[j]
        bw_matrix[l, 6] = kernz[k]
        bw_matrix[l, 7] = (10:100)[l]
        
      }
      
      kerns_list[[k]] = bw_matrix %>% 
        as.data.frame %>% 
        `colnames<-` (c("est", "se", "pv", "outcomes", "poly", "kern", "bw")) %>% 
        mutate(est = as.numeric(as.character(est)),
               se = as.numeric(as.character(se)),
               pv = as.numeric(as.character(pv)),
               outcomes = as.character(outcomes),
               poly = as.numeric(as.character(poly)),
               kern = as.character(kern),
               bw = as.numeric(as.character(bw)))
      
    }
    
    polys_list[[j]] = kerns_list
      
  }
  
  dataset_list[[i]] = polys_list
  
}

cct_estimates = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(kernpoly = paste0(kern, ", ", poly)) %>% 
  mutate(kernpoly = factor(kernpoly, 
                           levels = c("tri, 0", "tri, 1", "tri, 2", "tri, 3",
                                      "uni, 0", "uni, 1", "uni, 2", "uni, 3",
                                      "epa, 0", "epa, 1", "epa, 2", "epa, 3"),
                           labels = c("Tri., 0", "Tri., 1", "Tri., 2", "Tri., 3",
                                      "Uni., 0", "Uni., 1", "Uni., 2", "Uni., 3",
                                      "Epa., 0", "Epa., 1", "Epa., 2", "Epa., 3"))) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_stops", "count_calls", "hit", "arrest",
                                      "response_time", "rr_bw", "cat_person",
                                      "cat_property", "cat_society"),
                           labels = c("A. Terry Stops", "B. Officer-Initiated Calls", "C. Hit Rates",
                                      "D. Arrest Rates", "E. Response Time", "F. B/W Rate Ratio",
                                      "G. Crime (Against Person)", "H. Crime (Property)", "I. Crime (Against Society)"))) 

the_outs = c("A. Terry Stops", "B. Officer-Initiated Calls", "C. Hit Rates",
  "D. Arrest Rates", "E. Response Time", "F. B/W Rate Ratio",
  "G. Crime (Against Person)", "H. Crime (Property)", "I. Crime (Against Society)")

the_outs_labs = c("Terry Stops", "Officer-Initiated Calls", "Hit Rates",
             "Arrest Rates", "Response Time", "B/W Rate Ratio",
             "Crime (Against Person)", "Crime (Property)", "Crime (Against Society)")

the_outs_plot_list = as.list(rep(NA, length(the_outs_labs)))

for (i in 1:length(the_outs_labs)) {
  
  print(paste0("Iteration ", i))
  
  the_outs_plot_list[[i]] = 
    cct_estimates %>%
    filter(outcomes == the_outs[i]) %>% 
    ggplot() +
    geom_point(aes(x = bw, y = est),
               size = .1,
               alpha = .4) +
    geom_errorbar(aes(x = bw, 
                      ymin = est - 1.96 * se,
                      ymax = est + 1.96 * se),
                  width = 0,
                  size = .1) +
    geom_errorbar(aes(x = bw, 
                      ymin = est - 1.645 * se,
                      ymax = est + 1.645 * se),
                  width = 0,
                  size = .2) + 
    geom_hline(yintercept = 0) + 
    facet_wrap(~ kernpoly, scales = "free") + 
    labs(x = "Bandwidth", y = "RDiT Coefficient\n(BLM Protest)",
         title = the_outs_labs[i]) + 
    theme_tufte(base_size = 9)
  
}

ggsave(plot = the_outs_plot_list[[1]], filename = "pics/c2bwp1.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[2]], filename = "pics/c2bwp2.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[3]], filename = "pics/c2bwp3.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[4]], filename = "pics/c2bwp4.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[5]], filename = "pics/c2bwp5.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[6]], filename = "pics/c2bwp6.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[7]], filename = "pics/c2bwp7.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[8]], filename = "pics/c2bwp8.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[9]], filename = "pics/c2bwp9.png", width = 8, height = 6.25)

#### rdit: shaving 100 days after the protest #### 

# loop 

polys = c(0, 1, 2, 3)
kernz = c("tri", "uni", "epa")
datasets = list(terry_cct %>% as.data.frame,
                call_cct %>% as.data.frame,
                terry_hit_cct %>% as.data.frame,
                terry_arr_cct %>% as.data.frame, 
                call_rt_cct %>% as.data.frame,
                terry_bw_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame)
outcomes = 
  c("count_stops", "count_calls", "hit", 'arrest', 'response_time',
    "rr_bw", 'cat_person', "cat_property", "cat_society")

dataset_list = as.list(rep(NA, length(datasets)))

for (i in 1:length(datasets)) {
  
  print(paste0("I iteration ", i))
  polys_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J iteration ", j))
    
    # kerns_matrix = matrix(NA, ncol = 6, nrow = length(kernz))
    kerns_list = as.list(rep(NA, length(kernz)))
    
    for (k in 1:length(kernz)) {
      
      print(paste0("K iteration ", k))
      
      bw_matrix = matrix(NA, ncol = 7, nrow = length(1:100))
      
      for (l in 1:length(1:100)) {
        
        print(paste0("L iteration ", l))
        dftouse = 
          datasets[[i]] %>% 
          filter(blmrv <= -1 | blmrv >= l)
        
        dftouse$blmrv = 
          ifelse(dftouse$blmrv >= 0, dftouse$blmrv - l, dftouse$blmrv)
        
        out = rdrobust(y = as.numeric(dftouse[, outcomes[i]]),
                       x = dftouse$blmrv, 
                       p = polys[j],
                       kernel = kernz[k])
        
        bw_matrix[l, 1] = out$coef[3]
        bw_matrix[l, 2] = out$se[3]
        bw_matrix[l, 3] = out$pv[3]
        bw_matrix[l, 4] = outcomes[i]
        bw_matrix[l, 5] = polys[j]
        bw_matrix[l, 6] = kernz[k]
        bw_matrix[l, 7] = (1:100)[l]
        
      }
      
      kerns_list[[k]] = bw_matrix %>% 
        as.data.frame %>% 
        `colnames<-` (c("est", "se", "pv", "outcomes", "poly", "kern", "days_cut")) %>% 
        mutate(est = as.numeric(as.character(est)),
               se = as.numeric(as.character(se)),
               pv = as.numeric(as.character(pv)),
               outcomes = as.character(outcomes),
               poly = as.numeric(as.character(poly)),
               kern = as.character(kern),
               days_cut = as.numeric(as.character(days_cut)))
      
    }
    
    polys_list[[j]] = kerns_list
    
  }
  
  dataset_list[[i]] = polys_list
  
}

cct_estimates = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(kernpoly = paste0(kern, ", ", poly)) %>% 
  mutate(kernpoly = factor(kernpoly, 
                           levels = c("tri, 0", "tri, 1", "tri, 2", "tri, 3",
                                      "uni, 0", "uni, 1", "uni, 2", "uni, 3",
                                      "epa, 0", "epa, 1", "epa, 2", "epa, 3"),
                           labels = c("Tri., 0", "Tri., 1", "Tri., 2", "Tri., 3",
                                      "Uni., 0", "Uni., 1", "Uni., 2", "Uni., 3",
                                      "Epa., 0", "Epa., 1", "Epa., 2", "Epa., 3"))) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_stops", "count_calls", "hit", "arrest",
                                      "response_time", "rr_bw", "cat_person",
                                      "cat_property", "cat_society"),
                           labels = c("A. Terry Stops", "B. Officer-Initiated Calls", "C. Hit Rates",
                                      "D. Arrest Rates", "E. Response Time", "F. B/W Rate Ratio",
                                      "G. Crime (Against Person)", "H. Crime (Property)", "I. Crime (Against Society)"))) 

the_outs = c("A. Terry Stops", "B. Officer-Initiated Calls", "C. Hit Rates",
             "D. Arrest Rates", "E. Response Time", "F. B/W Rate Ratio",
             "G. Crime (Against Person)", "H. Crime (Property)", "I. Crime (Against Society)")

the_outs_labs = c("Terry Stops", "Officer-Initiated Calls", "Hit Rates",
                  "Arrest Rates", "Response Time", "B/W Rate Ratio",
                  "Crime (Against Person)", "Crime (Property)", "Crime (Against Society)")

the_outs_plot_list = as.list(rep(NA, length(the_outs_labs)))

for (i in 1:length(the_outs_labs)) {
  
  print(paste0("Iteration ", i))
  
  the_outs_plot_list[[i]] = 
    cct_estimates %>%
    filter(outcomes == the_outs[i]) %>% 
    ggplot() +
    geom_point(aes(x = days_cut, y = est),
               size = .1,
               alpha = .4) +
    geom_errorbar(aes(x = days_cut, 
                      ymin = est - 1.96 * se,
                      ymax = est + 1.96 * se),
                  width = 0,
                  size = .1) +
    geom_errorbar(aes(x = days_cut, 
                      ymin = est - 1.645 * se,
                      ymax = est + 1.645 * se),
                  width = 0,
                  size = .2) + 
    geom_hline(yintercept = 0) + 
    facet_wrap(~ kernpoly, scales = "free") + 
    labs(x = "Days Cut", y = "RDiT Coefficient\n(BLM Protest)",
         title = the_outs_labs[i]) + 
    theme_tufte(base_size = 9)
  
}

ggsave(plot = the_outs_plot_list[[1]], filename = "pics/dcp1.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[2]], filename = "pics/dcp2.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[3]], filename = "pics/dcp3.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[4]], filename = "pics/dcp4.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[5]], filename = "pics/dcp5.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[6]], filename = "pics/dcp6.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[7]], filename = "pics/dcp7.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[8]], filename = "pics/dcp8.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[9]], filename = "pics/dcp9.png", width = 8, height = 6.25)

#### temporal placebo tests #### 

# setting up loop

polys = c(0, 1, 2, 3)
kernz = c("uni")
datasets = list(terry_cct %>% as.data.frame,
               call_cct %>% as.data.frame,
               terry_hit_cct %>% as.data.frame,
               terry_arr_cct %>% as.data.frame,
               call_rt_cct %>% as.data.frame,
               terry_bw_cct %>% as.data.frame,
               crime_cct %>% as.data.frame,
               crime_cct %>% as.data.frame,
               crime_cct %>% as.data.frame)

outcomes =
 c("count_stops", "count_calls", "hit", 'arrest', 'response_time',
   "rr_bw", 'cat_person', "cat_property", "cat_society")

dataset_list = as.list(rep(NA, length(datasets)))

the_lag_max = as.numeric(as.Date("2020-05-29") - min(datasets[[1]]$date)) - 30

# loop

for (i in 1:length(datasets)) {

 print(paste0("I iteration ", i))
 polys_list = as.list(rep(NA, length(polys)))

 for (j in 1:length(polys)) {

   print(paste0("J iteration ", j))

   kerns_list = as.list(rep(NA, length(kernz)))

   for (k in 1:length(kernz)) {

     print(paste0("K iteration ", j))
     the_lag_max = (as.numeric(as.Date("2020-05-29") -
                                 min(datasets[[i]]$date))) - 30
     the_lag_max2 = as.Date(min(datasets[[i]]$date) + 30)
     lagmat = matrix(NA, ncol = 9, nrow = the_lag_max) %>%
       as.data.frame %>%
       `colnames<-` (c("est", "se", "pv", "outcomes", "poly",
                       "kern", "lag", "lagmax", "lagmaxdate"))

     for (z in 30:the_lag_max) {

       print(paste0("Z iteration ", z))

       datasets[[i]]$blmrv_fake = datasets[[i]]$date -
         datasets[[i]]$date[lead(datasets[[i]]$blmrv, z) == 0 &
                              !is.na(lead(datasets[[i]]$blmrv, z))]

       out = rdrobust(y = as.numeric(datasets[[i]][, outcomes[i]]),
                      x = datasets[[i]]$blmrv_fake,
                      p = polys[j],
                      kernel = kernz[k])

       lagmat[z, 1] = out$coef[3]
       lagmat[z, 2] = out$se[3]
       lagmat[z, 3] = out$pv[3]
       lagmat[z, 4] = outcomes[i]
       lagmat[z, 5] = polys[j]
       lagmat[z, 6] = kernz[k]
       lagmat[z, 7] = z
       lagmat[z, 8] = the_lag_max
       lagmat$lagmaxdate[z] = as.Date(the_lag_max2)

     }

     kerns_list[[k]] = lagmat %>% filter(!is.na(est)) %>%
       mutate(lagmaxdate = as.Date(lagmaxdate, origin = "1970-01-01"))

   }

   polys_list[[j]] =
     kerns_list

 }

 dataset_list[[i]] = polys_list

}

plc_test_out_sea = dataset_list
save(x = plc_test_out_sea, file = "data/plctest_sea.RData")
load(file = "data/plctest_sea.RData")

# re-running loop without placebos to serve as benchmarks

# loop 

polys = c(0, 1, 2, 3)
kernz = c("tri", "uni", "epa")
datasets = list(terry_cct %>% as.data.frame,
                call_cct %>% as.data.frame,
                terry_hit_cct %>% as.data.frame,
                terry_arr_cct %>% as.data.frame, 
                call_rt_cct %>% as.data.frame,
                terry_bw_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame)
outcomes = 
  c("count_stops", "count_calls", "hit", 'arrest', 'response_time',
    "rr_bw", 'cat_person', "cat_property", "cat_society")

dataset_list = as.list(rep(NA, length(datasets)))

for (i in 1:length(datasets)) {
  
  print(paste0("I iteration ", i))
  polys_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J iteration ", j))
    
    kerns_matrix = matrix(NA, ncol = 6, nrow = length(kernz))
    
    for (k in 1:length(kernz)) {
      
      print(paste0("K iteration ", k))
      out = rdrobust(y = as.numeric(datasets[[i]][, outcomes[i]]),
                     x = datasets[[i]]$blmrv, 
                     p = polys[j],
                     kernel = kernz[k])
      
      kerns_matrix[k, 1] = out$coef[3]
      kerns_matrix[k, 2] = out$se[3]
      kerns_matrix[k, 3] = out$pv[3]
      kerns_matrix[k, 4] = outcomes[i]
      kerns_matrix[k, 5] = polys[j]
      kerns_matrix[k, 6] = kernz[k]
      
    }
    
    polys_list[[j]] = 
      kerns_matrix %>% 
      as.data.frame %>% 
      `colnames<-` (c("est", "se", "pv", "outcomes", "poly", "kern")) %>% 
      mutate(est = as.numeric(as.character(est)),
             se = as.numeric(as.character(se)),
             pv = as.numeric(as.character(pv)),
             outcomes = as.character(outcomes),
             poly = as.numeric(as.character(poly)),
             kern = as.character(kern))
    
  }
  
  dataset_list[[i]] = polys_list
  
}

cct_estimates = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(kernpoly = paste0(kern, ", ", poly)) %>% 
  mutate(kernpoly = factor(kernpoly, 
                           levels = c("tri, 0", "tri, 1", "tri, 2", "tri, 3",
                                      "uni, 0", "uni, 1", "uni, 2", "uni, 3",
                                      "epa, 0", "epa, 1", "epa, 2", "epa, 3"),
                           labels = c("Tri.,\n0", "Tri.,\n1", "Tri.,\n2", "Tri.,\n3",
                                      "Uni.,\n0", "Uni.,\n1", "Uni.,\n2", "Uni.,\n3",
                                      "Epa.,\n0", "Epa.,\n1", "Epa.,\n2", "Epa.,\n3"))) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_stops", "count_calls", "hit", "arrest",
                                      "response_time", "rr_bw", "cat_person",
                                      "cat_property", "cat_society"),
                           labels = c("A. Terry Stops", "B. Officer-Initiated Calls", "C. Hit Rates",
                                      "D. Arrest Rates", "E. Response Time", "F. B/W Rate Ratio",
                                      "G. Crime (Against Person)", "H. Crime (Property)",
                                      "I. Crime (Against Society)"))) 

# placebo test plots 

plc_test_out_sea_df = bind_rows(
  plc_test_out_sea[[1]] %>% bind_rows(),
  plc_test_out_sea[[2]] %>% bind_rows(),
  plc_test_out_sea[[3]] %>% bind_rows(),
  plc_test_out_sea[[4]] %>% bind_rows(),
  plc_test_out_sea[[5]] %>% bind_rows(),
  plc_test_out_sea[[6]] %>% bind_rows(),
  plc_test_out_sea[[7]] %>% bind_rows(),
  plc_test_out_sea[[8]] %>% bind_rows(),
  plc_test_out_sea[[9]] %>% bind_rows()
) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_stops", "count_calls", "hit", 
                                      "arrest", "response_time", "rr_bw",
                                      "cat_person", "cat_property", "cat_society"),
                           labels = c("Terry Stops", "Officer Calls", "Hit Rate",
                                      "Arrest Rate", "Response Time", 
                                      "B/W Rate Ratio (Terry)", "Violent Crime",
                                      "Property Crime", "Society Crime")))

plc_test_out_sea_df$outcomes_poly = 
  factor(paste0(plc_test_out_sea_df$outcomes, ", ", plc_test_out_sea_df$poly),
         levels = unique(paste0(plc_test_out_sea_df$outcomes, ", ",
                                 plc_test_out_sea_df$poly)))

plc_test_out_sea_df2 = plc_test_out_sea_df %>% 
  filter(outcomes != "Response Time")
plc_test_out_sea_df2$outcomes_poly = 
  factor(plc_test_out_sea_df2$outcomes_poly, 
       levels = unique(plc_test_out_sea_df2$outcomes_poly))

cct_estimates2 = 
  (cct_estimates[grepl(x = cct_estimates$kernpoly, "Uni."), ] %>% 
  filter(outcomes != "E. Response Time"))

cct_estimates2$outcomes_poly = unique(plc_test_out_sea_df2$outcomes_poly)

# before plotting, figuring out the p-value

outcomes_poly_plc = 
  rep(NA, length(unique(plc_test_out_sea_df2$outcomes_poly)))

for (i in 1:length(unique(plc_test_out_sea_df2$outcomes_poly))) {
  
  print(paste0("Iteration ", i))
  
  outcomes_poly_plc[i] = 
    mean(abs(cct_estimates2$est[cct_estimates2$outcomes_poly == 
                                  unique(cct_estimates2$outcomes_poly)[i]]) >= 
           abs(plc_test_out_sea_df2$est[plc_test_out_sea_df2$outcomes_poly == 
                                          unique(plc_test_out_sea_df2$outcomes_poly)[i]]))
  
}

plc_test_out_sea_df2$outcomes_poly = 
  factor(plc_test_out_sea_df2$outcomes_poly,
         levels = unique(plc_test_out_sea_df2$outcomes_poly),
         labels = paste0(unique(plc_test_out_sea_df2$outcomes_poly), ", ",
                         round(outcomes_poly_plc, 2)))

cct_estimates2$outcomes_poly = 
  factor(cct_estimates2$outcomes_poly,
       levels = unique(cct_estimates2$outcomes_poly),
       labels = paste0(unique(cct_estimates2$outcomes_poly), ", ",
                       round(outcomes_poly_plc, 2)))

plc_test_out_sea_df2 %>% 
  ggplot() + 
  geom_histogram(aes(x = est)) + 
  facet_wrap(~ outcomes_poly, ncol = 4, scales = "free") +
  geom_vline(data = cct_estimates2,
             mapping = aes(xintercept = est)) +
  labs(x = "Placebo Coefficients", y = "Count") +  
  theme_tufte(base_size = 8)

ggsave(plot = last_plot(),
       filename = "pics/temp_plc_test_sea.png",
       width = 8, height = 8)

#### saving all datasets ####

datasets_sea = list(terry_cct %>% as.data.frame,
                call_cct %>% as.data.frame,
                terry_hit_cct %>% as.data.frame,
                terry_arr_cct %>% as.data.frame, 
                call_rt_cct %>% as.data.frame,
                terry_bw_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame,
                crime_cct %>% as.data.frame)

save(x = datasets_sea, file = 'data/datasets_sea.Rdata')
